#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "val_est.h"
#include "myrandom.h"
#include "func.h"
#include "est.h"
#include "va.h"

static double T(double *N, unsigned int m, double *p, unsigned int n);
static double *create_N(double *ma, unsigned int m, double a, double b, unsigned int start_inter, unsigned int inter_quant);
static double *create_P(Fun_Distro f, double a, double b, unsigned int start_inter, unsigned int inter_quant);
static double *create_ma(unsigned int id, unsigned int n, double param1, double param2);
static void calc_estimators(unsigned int id, double *ma_aux, unsigned int n);
static double g1Sim(double argv);
static double g2Sim(double argv);
static double g3Sim(double argv);


double Param1Sim = 0.0; /*valores globales para los parametros simulados de las tres distribuciones*/ 
double Param2Sim = 0.0;

double valor_p_chi_square(Fun_Distro f, double *ma, unsigned int n, unsigned int inter_quant, double a, double b, unsigned int start){
	double t=0.0, valor_p = 0.0;
	double *N = NULL, *P=NULL;
	
	N = create_N(ma,n,a,b,start,inter_quant);
	P = create_P(f,a,b,start,inter_quant);
	t = T(N,inter_quant,P,n);

	printf("\tvalor de t = %f\n",t);
	valor_p = chi_squareq(inter_quant-1.0,t);
	
	free(N);
	free(P);
	
	return valor_p;
}

/*estima el valor_p usando uniformes*/ 
double valor_p_KS_u(unsigned int n, double d0, unsigned int Num_Iter){
	unsigned int i = 0, j = 0, sum = 0;
	double *ma = NULL, d = 0.0;
	
	ma = (double *)calloc(n, sizeof(double));
	
	for(i = 0 ; i < Num_Iter ; i++){
		for(j = 0 ; j < n; j++)
			ma[j] = (1.0/(float)ULONG_MAX) * mzran13();
		
		bizzard_sort(ma, n);
		
		d = D(U1, ma, n);
		
		if(d > d0)
			sum++;
	}
	
	free(ma);
	ma = NULL;
	
	return (double)sum / (double) Num_Iter;
}

/*estima el valor_p simulando. Sirve para las 3 distribuciones. Id sirve para identificar de cual se trata*/
double valor_p_KS_sim(double d0,unsigned int id, double *ma, unsigned int n, double param1,double param2, unsigned int Num_Iter){
	unsigned int i = 0;
	double *ma_aux = NULL;
	double d=0.0, sum=0.0;
	
  
	for(i = 0; i<Num_Iter; i++){
		ma_aux = create_ma(id,n,param1, param2);/*creo muestra simulada con los parametros estimados de la muestra original*/
		calc_estimators(id,ma_aux,n); /*Seteo las variables globales Param1Sim y Param2Sim segun la muestra simulada*/
		bizzard_sort(ma_aux,n);
		
		/*valor de d simulado*/
		if(id == 0)
			d = D(g1Sim,ma_aux,n);
		else if(id == 1)	
			d = D(g2Sim,ma_aux,n);
		else
			d = D(g3Sim,ma_aux,n);
		 
		if(d >= d0)
			sum++;
		  
		free(ma_aux);
	}
	
	return sum/(double)Num_Iter;
}

double D(Acu_Fun f, double *ma, unsigned int n){
	unsigned int i = 0;
	double res = 0.0, der = 0.0, izq = 0.0;
	
	for(i = 1 ; i < n ; i++){
		double F = f(ma[i-1]);
		
		izq = ((double)i / (double)n) - F;
		der = F - (((double)i - 1) / (double)n);
		
		if(izq > der && izq > res)
			res = izq;
		else if(der > izq && der > res)
			res = der;
	}
	return res;
}

/*Crea muestra simulada para las 3 distribuciones*/
static double *create_ma(unsigned int id, unsigned int n, double param1, double param2){
	unsigned int i = 0;
	double *ma = calloc(n,sizeof(double));
	double *res = NULL;
	
	if(id == 0){ /*genero muestra aleatoria normal*/
		while(i < n){
			res = Normal_MP(param1, param2); /*genera 2 normales*/
			ma[i] = res[0];
			i++;
			if( i < n)
				ma[i] = res[1];
			i++;
			free(res);
		}
	}
	else if(id ==1){/*genero muestra aleatoria lognormal*/
		for(i = 0; i < n; i++)
			ma[i]= lognormal(param1, param2);
	}
	else{ /*genero muestra aleatoria gamma*/
		for(i = 0; i < n; i++)
			ma[i] = gamma(param1, param2);
	}
			
	return ma;
}

/*Setea los variables globales*/
static void calc_estimators(unsigned int id, double *ma_aux, unsigned int n){
	
	if(id == 0){/*estimo mu y sigma para la normal*/
		Param1Sim = Estimador_MU_Normal(ma_aux,n);
		Param2Sim = pow(Estimador_V_Normal(ma_aux,n),0.5);
	}
	else if(id == 1){/*estimo mu y varianza para la lognormal*/
		Param1Sim = Estimador_MU_Log_Normal(ma_aux, n);
		Param2Sim = pow(Estimador_V_Log_Normal(ma_aux,n),0.5);
	}
	else{/*estimo alpha y lambda para la gamma**/
		Param1Sim = Estimador_Alpha_Gamma(ma_aux,n);
		Param2Sim = Estimador_Lambda_Gamma(Param1Sim, ma_aux,n);
	}
}

/*Son las F's con los parametros simulados*/
static double g1Sim(double argv){
	return normalp(argv, Param1Sim, Param2Sim);
}

static double g2Sim(double argv){
	return log_normalp(argv, Param1Sim, Param2Sim);
}

static double g3Sim(double argv){
	return gammap(argv, Param1Sim, 1.0/Param2Sim);
}

static double T(double *N, unsigned int m, double *p, unsigned int n){
        unsigned int i = 0;
        double T = 0.0;
	double np = 1.0;
	double temp = 0.0;
        
	for(i=0; i < m; i++){
		if(p[i]!=0.0){
			np = (double)n*p[i];
			temp = pow(N[i] - np, 2.0)/np;
			T = T + temp;
		}
	}
	 
        return T;
}

static double *create_N(double *ma, unsigned int m, double a, double b, unsigned int start_inter, unsigned int inter_quant){
        unsigned int i = 0, j = 0;
        double *N = NULL;
        double inter_size = (a + b)/inter_quant;
        
        N = (double *)calloc(inter_quant, sizeof(double));
        
        
        for(i = 0 ; i < inter_quant ; i++){
                for(j = 0 ; j < m ; j++){
                        if((start_inter + ((double)i)*inter_size) <= ma[j] &&  ma[j] < start_inter + ((double)i + 1.0)*inter_size)
                                N[i]++;
                }
        }
        return N;
}

static double *create_P(Fun_Distro f, double a, double b, unsigned int start_inter, unsigned int inter_quant){
        double *P = NULL, inter_size = (a + b)/inter_quant, sum = 0.0, j = 0.0, sum2 = 0.0;
        unsigned int i = 0;
        
        P = (double *)calloc(inter_quant, sizeof(double));
        	
	for(i = 0 ; i < inter_quant-1 ; i++){
	    for(j = start_inter + (double)i*inter_size; j <= (double)(start_inter + i)*inter_size;j++)
		sum+=f(j);
	    
	    P[i] = sum;
	    sum2+=sum;
	    sum = 0;
	}
	
	P[i] = 1.0 - sum2;
       
        return P;
}

